#!/usr/bin/env julia

"""
Basic 2D Lid-Driven Cavity Example

This example demonstrates the basic usage of LidDrivenCavity.jl for solving
a 2D lid-driven cavity flow on an L-shaped domain.

Run with: julia basic_2d_example.jl
"""

# Add the package to the path
push!(LOAD_PATH, "../src")

using LidDrivenCavity
using Printf

function main()
    println("🌀 LidDrivenCavity.jl - Basic 2D Example")
    println("="^50)
    
    # Display package information
    lid_driven_info()
    println()
    
    # Create solver options
    options = LidDrivenOptions(
        n = 16,                    # Grid size (16x16 interior points)
        Re = 100.0,               # Reynolds number
        max_steps = 1000,         # Maximum iterations
        tol = 1e-5,               # Convergence tolerance
        solver = :julia,          # Use Julia backend 
        preconditioner = :diagonal,
        verbose = true,           # Detailed output
        save_history = true       # Save convergence history
    )
    
    println("🔧 Solver Configuration:")
    println("  Grid size: $(options.n) × $(options.n)")
    println("  Reynolds number: $(options.Re)")
    println("  Max iterations: $(options.max_steps)")
    println("  Tolerance: $(options.tol)")
    println("  Solver: $(options.solver)")
    println()
    
    # Solve the problem
    println("🚀 Starting solver...")
    t_start = time()
    
    try
        result = solve_lid_driven_2d(options)
        
        solve_wall_time = time() - t_start
        
        println("\n✅ Solution completed!")
        println("="^50)
        
        # Print results summary
        println("📊 Results Summary:")
        println("  Converged: $(result.converged)")
        println("  Iterations: $(result.iterations)")
        println("  Final residual: $(@sprintf("%.2e", result.residual_norm))")
        println("  Solver time: $(@sprintf("%.3f", result.solve_time)) seconds")
        println("  Wall time: $(@sprintf("%.3f", solve_wall_time)) seconds")
        
        if !isempty(result.convergence_history)
            println("  Initial residual: $(@sprintf("%.2e", result.convergence_history[1]))")
            reduction = result.convergence_history[1] / result.residual_norm
            println("  Residual reduction: $(@sprintf("%.1e", reduction))")
        end
        
        # Solution statistics
        println("\n🌊 Flow Field Statistics:")
        ψ_range = extrema(result.ψ[.!isnan.(result.ψ)])
        ω_range = extrema(result.ω[.!isnan.(result.ω)])
        u_max = maximum(abs.(result.u))
        v_max = maximum(abs.(result.v))
        
        println("  Stream function ψ: [$(@sprintf("%.3f", ψ_range[1])), $(@sprintf("%.3f", ψ_range[2]))]")
        println("  Vorticity ω: [$(@sprintf("%.3f", ω_range[1])), $(@sprintf("%.3f", ω_range[2]))]")
        println("  Max |u|: $(@sprintf("%.3f", u_max))")
        println("  Max |v|: $(@sprintf("%.3f", v_max))")
        
        # Convergence analysis
        if !isempty(result.convergence_history)
            conv_analysis = convergence_analysis(result)
            if conv_analysis.analysis_available
                println("\n📈 Convergence Analysis:")
                println("  Convergence rate: $(@sprintf("%.3f", conv_analysis.convergence_rate))")
                println("  Reduction factor: $(@sprintf("%.1e", conv_analysis.reduction_factor))")
            end
        end
        
        println("\n💾 Solution data available in 'result' object")
        println("   Fields: ψ (stream function), ω (vorticity), u, v (velocities)")
        println("   Grid: x, y coordinates")
        println("   Options: solver configuration used")
        
        return result
        
    catch e
        println("\n❌ Solver failed with error:")
        println("   $e")
        
        if isa(e, LoadError) || isa(e, MethodError)
            println("\n💡 Possible solutions:")
            println("   • Make sure all dependencies are available")
            println("   • Try with smaller grid size: n=8")  
            println("   • Check that GCR.jl is properly installed")
        end
        
        return nothing
    end
end

# Run the example
if abspath(PROGRAM_FILE) == @__FILE__
    result = main()
    
    if result !== nothing
        println("\n🎉 Example completed successfully!")
        println("   Use plot_solution(result) to visualize (requires Plots.jl)")
    else
        println("\n😞 Example failed - check error messages above")
    end
end